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1  Summary 

During  the  past  two  years  the  investigators  have  been  able  to  develop  a  computer  code  which, 
on  the  basis  of  some  preliminary  experiments,  has  turned  out  to  be  quite  competitive  with  a  well 
established  code.  The  algorithmic  development  phase  of  the  research  has  led  to  several  useful 
results.  An  eighth  order  method,  which  has  second,  fourth  and  sixth  order  methods  embedded 
in  it  has  been  developed.  The  method  possesses  some  of  the  best  features  of  implicit  Runge- 
Kutta  and  Gap  schemes.  A  multiderivative  generalization  of  the  above  schemes  has  also  been 
realized.  A  sparse  factorization  for  quasi-Newton  type  methods  has  been  obtained.  Algorithms 
for  solving  combined  systems  of  linear  and  non-linear  algebraic  equations  and  matrix  splitting 
have  also  been  developed.  A  new  algorithm  has  been  found  for  updating  the  sparse  LDU 
factorization  of  the  approximation  to  the  Jacobian  of  the  system  of  equations.  This  method 
is  relatively  more  efficient  than  the  currently  available  methods.  A  practical  updating  method 
has  been  obtained  that  is  in  the  Broyden  family  of  updates  and  is  in  most  cases  better  than 
the  BFGS  method.  Some  effort  was  -directed  towards  improvement  in  the  stability  properties 
of  Boundary  Value  Runge  Kuttar(BVRK^methods. 

Some  of  our  results  are  being  used  for  the  computer  simulation  of  kidney  function  at  the 
Cornell  Medical  School  in  a  collaborative  research  in  which  the  Principal  Investigator  is  involved. 


2  Report 

In  [lj,  we  consider  a  system  of  differential  equations 

y'(x)  =  f(x,y)  ,  x  e  [a, 6] 
with  two-point  boundary  conditions 


g{y{a),y(b))  =  0. 

For  the  subinterval  [x0,X|],  the  proposed  numerical  methods  can  be  written  as 


y,  =  (1  -  ©i)yo  +  0.yi  +  h  ^a0/(x0  +  Cjh,yj) 

j= i 


i  —  1, s 

3 

y,  =  y0  +  btf(x0  +  c,h,yj) 


i=i 


with 


at  J  —  0  V  j  >  i. 

We  will  call  them  Boundary  Value  Runge  Kutta  (BVRK)  processes.  These  finite  difference 
schemes  combine  features  of  both  Runge- Kutta  processes  and  Gap  schemes.  Just  like  Runge- 
Kutta  processes,  these  methods  do  not  require  derivatives  of  function  /.  Also,  just  like  Gap 


U 

□ 


ides 

or 


schemes,  BVRK  methods  use  information  available  at  both  end-points  of  a  subinterval.  No 
information  available  outside  a  subinterval  is  used  for  approximating  the  ODE  in  that  subin¬ 
terval.  Because  of  symmetry,  the  asymptotic  error  expansion  for  the  numerical  solution  has 
only  even  powers  of  mesh-size  h.  It  is  possible  to  develop  A-stabie  BVRK  processes,  which  is 
desirable.  We  have  developed  some  techniques  which  can  be  used  to  derive  BVRK  processes. 
A  quasi- Newton  method  for  the  solution  of  the  systems  of  nonlinear  equations  arising  from 
the  use  of  these  methods,  has  also  been  developed.  We  have  derived  an  eighth  order  BVRK 
method  which  has  second,  fourth,  and  sixth  embedded  in  it.  We  have  developed  a  technique 
which  takes  advantage  of  this  natural  embedding  of  methods  for  generating  asymptotically 
equidistributing  grids,  and  for  estimating  the  error  in  the  numerical  solution.  We  have  devel¬ 
oped  a  general  purpose  adaptive  solver  based  on  these  ideas.  In  the  future,  we  hope  to  make 
this  code  compatible  with  PASVA3  (which  is  a  popular  general  purpose  adaptive  code)  so  that 
depending  on  various  conditions  the  final  program  will  either  use  BVRK  processes  or  Deferred 
correction. 

In  (2],  we  consider  a  multiderivative  generalization  of  BVRK  processes.  These  methods 
use  fewer  derivatives  of  function  /  to  yield  results  which  are  the  same  as  the  results  for  Gap 
schemes.  Hence,  these  processes  will  be  preferred  over  Gap  schemes  if  the  higher  order  analytic 
derivatives  are  not  available. 

A  general  description  of  sparsity  based  algorithms  that  use  structural  information  in  handling 
the  model  equations  is  given  in  [3].  The  primary  focus  is  on  splitting  of  equations  and  variables 
into  two  subsets:  one  large  but  relatively  essy  to  solve,  called  the  non-basic  subset,  and  the 
other  small  but  difficult  to  handle,  called  the  basic  subset.  The  non-basic  subset  is  then  solved 
more  often  than  the  difficult  basic  set.  The  computer  storage  is  that  required  by  the  basic  set. 
Applications  in  flow  networks  and  energy-economy  models  are  briefly  described. 

A  new  quasi- Newton  method  for  sparse  matrices  arising  in  various  application  areas  is 
presented  in  [4].  In  this  work  sparse  triple  factorization  is  used  to  develop  a  rank  two  update. 

A  class  of  processes  [5]  have  been  found  for  solving  second  order  boundary  value  problems 
of  the  form  y”  =  f(x,y),  y(a)=A,  y(b)=B.  These  do  not  require  derivatives  of  f,  and  result 
in  nonlinear  equations  with  blocked  tridiagonal  Jacobians.  Specialized  techniques  have  been 
developed  to  derive  the  necessary  formulas.  The  computation  of  the  Jacobians  for  the  dis¬ 
cretized  systems  arising  from  the  use  of  these  formulas  is  expensive.  This  makes  the  ordinary 
Newton  iterations  impractical.  This  paper  contains  various  techniques  for  approximating  and 
splitting  the  Jacobians.  These  techniques  are  computationally  inexpensive  and  lead  to  excellent 
convergence. 

An  efficient  method  for  handling  systems  with  linear  and  non-linear  subsystems  is  given  in 
[6].  In  this  method  iterations  are  needed  only  for  the  solution  of  non-linear  subsystems. 

Cubic  and  quintic  splines  have  been  utilized  in  some  of  the  methods  for  solving  differential 
equations  developed  by  the  investigators  and  others  in  the  past  .  One  of  the  crucial  steps  in 
these  techniques  involves  the  solution  of  band  systems.  A  very  efficient  general  purpose  scheme 
for  handling  such  band  systems  [7j  has  been  developed.  This  scheme  is  also  very  useful  for 
block  banded  systems  arising  in  a  variety  of  applications. 

A  method  of  order  eight  that  uses  cubic  splines  on  quintic  splines  for  solving  differential 


equations  is  described  in  [8]. 

The  Principal  Investigator  has  given  invited  papers  based  on  some  of  the  research  resulting 
from  this  grant  at  two  international  conferences  (9,10).  In  [11],  we  consider  the  numerical 
solution  of  parabolic  partial  differential  equations  of  Schrodinger  type.  Most  conventional 
explicit  finite  difference  methods  are  unconditionally  unstable  for  these  problems.  Use  of  explicit 
Runge-Kutta  methods  with  method  of  lines  also  yields  unstable  results.  However,  use  of  an 
implicit  method  requires  solution  of  complex  systems  of  algebraic  equations.  Since,  very  little 
is  known  about  large  complex  systems  as  compared  to  real  systems,  explicit  schemes  would 
be  preferred  over  implicit  schemes.  Also,  most  explicit  schemes  can  easily  be  vectorized.  In 
[11],  we  consider  the  use  of  dissipative  terms  for  deriving  conditionally  stable  explicit  single 
step  schemes.  We  also  consider  use  of  a  pair  of  methods  in  an  alternating  fashion  for  getting 
more  accurate  results.  Most  of  the  methods  considered  in  this  paper  can  easily  be  extended 
for  more  general  Schrodinger  equations. 

In  [12],  we  consider  the  following  singularly  perturbed  two-point  boundary  value 

-eu"  +  a( x)  u  —  /(x,u) 

«(0)  =  go  ,  u(l)  =  gi- 

We  assume  that  a  and  /  are  sufficiently  smooth  functions  and  that  a(x)  >  a  >  0  on  [0,1]. 
Then  using  the  usual  splitting  one  can  prove  that  for  e  sufficiently  small,  the  solution  of  the 
above  problem  can  be  represented  as 

1  /•! 

u(x-,e)  =  A(x;e)  +  J3(x;e)exp[  /  a] 

t  J  X 

where  A  and  B,  and  their  derivatives  up  to  an  integer  K  can  be  bounded  on  [0,1]  independently 
of  e.  For  subintervals  [xi_i,Xi]  and  [xj,x,+i],  we  consider  finite  difference  schemes  of  the 
following  form: 

j 

(Jj  =  OjUi- 1  +  SjUi  +  PiU,+ 1  -  ^2ijkf{'Pk,Uk)  ,  j  =  1,—,  J 

k=  1 

J 

-I-  ao  Ui  4-  ct\Ui  +  l  —  T.  fiifj  *Pi ,  Uj ) 

}= 

where  coefficients  0j,  6Jt  p;,  7^,  f3j,  a_i,  a0,  and  07  are  chosen  so  that  the  scheme  is  exact 
on  a  class  of  functions 

1  1  rx 

1,  X  ,  . ,xn,  exp(- J  a),  ....,  xLexp(- J  a). 

First,  second  and  third  order  methods  of  this  type  have  been  tried  on  some  test  problems  and 
they  exhibit  the  same  uniform  order  of  convergence  as  predicted  by  the  theory. 

In  [13],  we  consider  linear  complementarity  problems  having  special  structures.  A  large 
number  of  free-boundary  problems  can  be  formulated  as  linear-complementary  problems  having 
the  properties  discussed  in  this  paper. 

In  [14],  a  new  definition  of  stability  that  characterizes  the  behavior  of  BVRK  methods  is 
being  developed. 
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